data{
  int N; 
  int K; 
  vector[N] y_obs;
  matrix[N,K] X;
  int N_o;
  int N_d;
  int<lower=1,upper=N_o> i_o[N];
  int<lower=1,upper=N_d> i_d[N];
  int N_year;
  int<lower=1,upper=N_year> i_year[N];
}

parameters{
  vector[K] beta; 
  real<lower=0> shape;
  vector<lower=0,upper=1>[N_year] pi;
  real<lower=0> theta0;
  vector[N_o] e_theta_0;
  real s_theta_0;
  }

transformed parameters{
  vector[N] scale;
  {
    real mean_e_theta_0;
    real sd_e_theta_0;
    matrix[N_o, N_d] theta;
    
    mean_e_theta_0 = mean(e_theta_0);
    sd_e_theta_0 = sd(e_theta_0);
    for(n_o in 1:N_o){
      for(n_d in 1:N_d){
        theta[n_o, n_d] = theta0 + (s_theta_0 * 
          ((e_theta_0[i_o[n_o]] - mean_e_theta_0)/sd_e_theta_0));
      }
    }
     for(n in 1:N){
        scale[n] = exp(X[n] * beta) * exp(theta[i_o[n], i_d[n]]); 
       }
    }
}
model{
  //shape ~ normal(0, 10);
  //beta ~ normal(0, 10);
  //pi ~ normal(0, 1);
  shape ~ cauchy(0, 2.5);
  beta ~ normal(0, 10);
  pi ~ beta(2, 2);
  theta0 ~ normal(0, 5);
  e_theta_0 ~ normal(0, 1);
  s_theta_0 ~ cauchy(0, 2.5);
  for(n in 1:N){
    if(y_obs[n] == 0){
    target+= bernoulli_lpmf(1 | pi[i_year[n]]);
  }
  else{
        target+= bernoulli_lpmf(0 | pi[i_year[n]]) + 
          weibull_lpdf(y_obs[n] | shape, scale[n]);  
      }

  }
}
generated quantities {
  vector[N] y_rep;
  vector[N] resids;
  {
   vector[N] y_pred;
  for(n in 1: N){
   y_pred[n] = X[n]*beta;
  }
  resids = exp(y_pred) - y_obs;
  for(n in 1:N){
    int zero_obs;
    zero_obs = bernoulli_rng(pi[i_year[n]]);
    if(zero_obs == 1){
      y_rep[n] = 0;
    }
    else{
     y_rep[n] = weibull_rng(shape, scale[n]);
    }
  }
  }
}



